* Dofile for generating estimates in figures for "Signaling Restraints"  
* codes tested with stata version 14 as well as version 16

clear all
set more off

* install "more_clarify" package to use "sumqoi" command: if you already have this, it's fine to skip this line.
ssc install more_clarify, replace

* set working directory to the main folder 
cd "[the main folder]"

*====================================*
* 1. Figure 2: Non-military Support from NSA (Transnational support)
*====================================*

* Load the data
use "commitment data master.dta", clear


* Run the main model (Model 2)
quietly melogit commitment_dum transconstsupp_dum strengthcent_num_4s  peacetalk_dum  rebstrength_num_5s terrcont_dum  ///
 rebextpart_dum secessionist  age   rebpresosts_dum i.rtypesup_num_3i i.rebpolwing_3i numopp dem_gwf time ///
 t t2 t3  || dyadid2:

	tab transconstsupp_dum if e(sample), matcell(obs) // extract the number of obs. for each category
	matrix list obs
	svmat obs, names(obs)

	* Predict Random effects "mu"
	
	tempfile temp
	predict ra if e(sample), reffects /* Random intercepts */

	gen keep = 1 if e(sample)
	drop if keep !=1

	save `temp'


	set seed 100
	mat b = e(b)
	mat V = e(V)

	drawnorm b_transconstsupp b_strengthcent b_peacetalk b_rebstrength b_terrcont   ///
	 b_rebextpart b_secessionist b_age b_rebpresosts b_rtypesup0 b_rtypesup1 b_rtypesup2 b_rebpolwing0 b_rebpolwing1 b_rebpolwing2 ///
	 b_numopp    b_dem  b_time b_t b_t2 b_t3 b_cons b_ra, mean(b) cov(V) n(1000) clear


	merge using `temp'

	
	set more off


	* Calculate the pred. prob. of commitment  across transnational support when the other covariables set to their observed values
	forvalues y = 0/1 {

	gen psup_`y'_mean = .
	forvalues i = 1/1000 {
	gen psup_`y'_`i' = invlogit((b_cons[`i']) + strengthcent_num_4s*b_strengthcent[`i'] + ///
	rebstrength_num_5s*b_rebstrength[`i'] + terrcont_dum*b_terrcont[`i'] + ///
	`y'*b_transconstsupp[`i'] + peacetalk_dum*b_peacetalk[`i'] +  ///
	rebextpart_dum*b_rebextpart[`i'] + secessionist*b_secessionist[`i'] +  age*b_age[`i'] + ///
	rebpresosts_dum*b_rebpresosts[`i'] + ///
	0.rtypesup_num_3i*b_rtypesup0[`i'] + 1.rtypesup_num_3i*b_rtypesup1[`i'] + 2.rtypesup_num_3i*b_rtypesup2[`i'] + ///
	0.rebpolwing_3i*b_rebpolwing0[`i'] + 1.rebpolwing_3i*b_rebpolwing1[`i'] + 2.rebpolwing_3i*b_rebpolwing2[`i'] + ///
	numopp*b_numopp[`i'] + 	dem_gwf*b_dem[`i'] + time*b_time[`i'] + t*b_t[`i'] + t2*b_t2[`i'] + t3*b_t3[`i'])

	summarize psup_`y'_`i', meanonly
	replace psup_`y'_mean = r(mean) in `i'

	drop psup_`y'_`i'
	}
	}


	* print all predicted probabilities 
	/* The command "sumqoi" nicely summarizes the quantity of interests.
	To use the command "sumqoi", you need to install "more clarify" program. 
	Please install the package as written above. */


	sumqoi psup_*


* generate a csv file with the estimates to make Figure 2 in R.
	return list
	
	matrix mean=r(mean)
	matrix lb95=r(lb)
	matrix ub95=r(ub)
	svmat mean, names(mean)
	svmat lb95, names(lb95)
	svmat ub95, names(ub95)

	foreach v in mean lb95 ub95 obs {
rename `v'1 `v'
}
	keep mean lb95 ub95 obs
		gen id = _n
		drop if mean==.
		outsheet id mean lb95 ub95 obs using "figures/psup for figure 2.csv" , comma nolabel replace	
		


*====================================*
* 2. Figure 3: Central Control Strength
*====================================*


* Load the data
use "commitment data master.dta", clear

* Run the main model (Model 2)
quietly melogit commitment_dum transconstsupp_dum strengthcent_num_4s  peacetalk_dum  rebstrength_num_5s terrcont_dum  ///
 rebextpart_dum secessionist  age   rebpresosts_dum i.rtypesup_num_3i i.rebpolwing_3i numopp dem_gwf time ///
 t t2 t3  || dyadid2:


	tab strengthcent_num_4s if e(sample), matcell(obs) // extract the number of obs. for each category
	matrix list obs
	svmat obs, names(obs)

	* Predict Random effects "mu"
	
	tempfile temp
	predict ra if e(sample), reffects /* Random intercepts */

	gen keep = 1 if e(sample)
	drop if keep !=1

	save `temp'


	set seed 100
	mat b = e(b)
	mat V = e(V)

	drawnorm b_transconstsupp b_strengthcent b_peacetalk b_rebstrength b_terrcont   ///
	 b_rebextpart b_secessionist b_age b_rebpresosts b_rtypesup0 b_rtypesup1 b_rtypesup2 b_rebpolwing0 b_rebpolwing1 b_rebpolwing2 ///
	 b_numopp    b_dem  b_time b_t b_t2 b_t3 b_cons b_ra, mean(b) cov(V) n(1000) clear


	merge using `temp'

	set more off
	
	* Calculate the pred. prob. of commitment across the levels of central control strength when the other covariables set to their observed values
	forvalues y = 0/3 {

	gen pccs_`y'_mean = .
	forvalues i = 1/1000 {
	gen pccs_`y'_`i' = invlogit((b_cons[`i']) + `y'*b_strengthcent[`i'] + ///
	rebstrength_num_5s*b_rebstrength[`i'] + terrcont_dum*b_terrcont[`i'] + ///
	transconstsupp_dum*b_transconstsupp[`i'] + peacetalk_dum*b_peacetalk[`i'] +  ///
	rebextpart_dum*b_rebextpart[`i'] + secessionist*b_secessionist[`i'] +  age*b_age[`i'] + ///
	rebpresosts_dum*b_rebpresosts[`i'] + ///
	0.rtypesup_num_3i*b_rtypesup0[`i'] + 1.rtypesup_num_3i*b_rtypesup1[`i'] + 2.rtypesup_num_3i*b_rtypesup2[`i'] + ///
	0.rebpolwing_3i*b_rebpolwing0[`i'] + 1.rebpolwing_3i*b_rebpolwing1[`i'] + 2.rebpolwing_3i*b_rebpolwing2[`i'] + ///
	numopp*b_numopp[`i'] + 	dem_gwf*b_dem[`i'] + time*b_time[`i'] + t*b_t[`i'] + t2*b_t2[`i'] + t3*b_t3[`i'])

	summarize pccs_`y'_`i', meanonly
	replace pccs_`y'_mean = r(mean) in `i'

	drop pccs_`y'_`i'
	}
	}


	* print all predicted probabilities 
	/* The command "sumqoi" nicely summarizes the quantity of interests.
	To use the command "sumqoi", you need to install "more clarify" program. 
	Please install the package as written above. */


* generate a csv file with the estimates to make Figure 3 in R.

	sumqoi pccs_*


	matrix mean=r(mean)
	matrix lb95=r(lb)
	matrix ub95=r(ub)
	svmat mean, names(mean)
	svmat lb95, names(lb95)
	svmat ub95, names(ub95)

	foreach v in mean lb95 ub95 obs {
rename `v'1 `v'
}
	keep mean lb95 ub95 obs
		gen id = _n
		drop if mean==.
		outsheet id mean lb95 ub95 obs using "figures/pccs for figure 3.csv" , comma nolabel replace	




*====================================*
* 3. Figure 4: Near peace talk 
*====================================*


* Load the data
use "commitment data master.dta", clear


* Run the main model (Model 2)
quietly melogit commitment_dum transconstsupp_dum strengthcent_num_4s  peacetalk_dum  rebstrength_num_5s terrcont_dum  ///
 rebextpart_dum secessionist  age   rebpresosts_dum i.rtypesup_num_3i i.rebpolwing_3i numopp dem_gwf time ///
 t t2 t3  || dyadid2:


	tab peacetalk_dum if e(sample), matcell(obs) // extract the number of obs. for each category
	matrix list obs
	svmat obs, names(obs)


	* Predict Random effects "mu"
	
	tempfile temp
	predict ra if e(sample), reffects /* Random intercepts */

	gen keep = 1 if e(sample)
	drop if keep !=1

	save `temp'


	set seed 100
	mat b = e(b)
	mat V = e(V)

	drawnorm b_transconstsupp b_strengthcent b_peacetalk b_rebstrength b_terrcont   ///
	 b_rebextpart b_secessionist b_age b_rebpresosts b_rtypesup0 b_rtypesup1 b_rtypesup2 b_rebpolwing0 b_rebpolwing1 b_rebpolwing2 ///
	 b_numopp    b_dem  b_time b_t b_t2 b_t3 b_cons b_ra, mean(b) cov(V) n(1000) clear


	merge using `temp'

	set more off

	* Calculate the pred. prob. of commitment  when the other covariables set to their observed values
	forvalues y = 0/1 {

	gen npt_`y'_mean = .
	forvalues i = 1/1000 {
	gen npt_`y'_`i' = invlogit((b_cons[`i']) + strengthcent_num_4s*b_strengthcent[`i'] + ///
	rebstrength_num_5s*b_rebstrength[`i'] + terrcont_dum*b_terrcont[`i'] + ///
	transconstsupp_dum*b_transconstsupp[`i'] + `y'*b_peacetalk[`i'] +  ///
	rebextpart_dum*b_rebextpart[`i'] + secessionist*b_secessionist[`i'] +  age*b_age[`i'] + ///
	rebpresosts_dum*b_rebpresosts[`i'] + ///
	0.rtypesup_num_3i*b_rtypesup0[`i'] + 1.rtypesup_num_3i*b_rtypesup1[`i'] + 2.rtypesup_num_3i*b_rtypesup2[`i'] + ///
	0.rebpolwing_3i*b_rebpolwing0[`i'] + 1.rebpolwing_3i*b_rebpolwing1[`i'] + 2.rebpolwing_3i*b_rebpolwing2[`i'] + ///
	numopp*b_numopp[`i'] + 	dem_gwf*b_dem[`i'] + time*b_time[`i'] + t*b_t[`i'] + t2*b_t2[`i'] + t3*b_t3[`i'])

	summarize npt_`y'_`i', meanonly
	replace npt_`y'_mean = r(mean) in `i'

	drop npt_`y'_`i'
	}
	}


	* print all predicted probabilities 
	/* The command "sumqoi" nicely summarizes the quantity of interests.
	To use the command "sumqoi", you need to install "more clarify" program. 
	Please install the package as written above. */

	sumqoi npt_*


* generate a csv file with the estimates to make Figure 4 in R.

	matrix mean=r(mean)
	matrix lb95=r(lb)
	matrix ub95=r(ub)
	svmat mean, names(mean)
	svmat lb95, names(lb95)
	svmat ub95, names(ub95)

	foreach v in mean lb95 ub95 obs {
rename `v'1 `v'
}
	keep mean lb95 ub95 obs
		gen id = _n
		drop if mean==.
		outsheet id mean lb95 ub95 obs  using "figures/pnpt for figure 4.csv" , comma nolabel replace	
